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ABSTRACT 

This report describes the physical models employed in GPSOMC, the modeling module of the 
GIPSY software system developed at JPL for analysis of geodetic Global Positioning Satellite (GPS) 
measurements. Details of the various contributions to range and phase observables are given, as well 
as the partial derivatives of the observed quantities with respect to model parameters. A glossary of 
parameters is provided to enable persons doing data analysis to identify quantities in the current report 
with their counterparts in the computer programs. The present version is a revision of the original 
document of the same title, JPL Publication 87-21, dated September 15, 1987, which it supersedes. 
There are no basic model revisions, with the exceptions of an improved ocean loading model, described 
in Section 2.2.2, and some new options for handling clock parametrization, mentioned in Section 3 
and Appendix B. Such misprints as were discovered have been corrected. The authors hope to publish 
further revisions of this document in the future, both to include modeling improvements, and to ensure 
that the model description is always in accord with the current software. 


iii 


CONTENTS 


1. INTRODUCTION 1 

2. COORDINATE FRAME AND GEOMETRY 2 

2.1 RECEIVER LOCATION TIME DEPENDENCE 2 

2.2 TIDAL EFFECTS 2 

2.2.1 Solid Earth Tides 3 

2.2.2 Ocean Loading 4 

2.2.3 Pole Tide 5 

2.3 TRANSFORMATION FROM EARTH-FIXED TO GEOCENTRIC INERTIAL 

COORDINATE SYSTEMS 5 

2.4 UT1 AND POLAR MOTION 6 

2.5 NUTATION 9 

2.6 PRECESSION 14 

2.7 PERTURBATION ROTATION 15 

2.8 GEOCENTER OFFSET AND COORDINATE SCALING 16 

2.9 PHASE CENTER OFFSETS 16 

3. OBSERVABLES AND CLOCK PARAMETERS 17 

3.1 PSEUDO-RANGE 18 

3.2 CARRIER PHASE 20 

4. TROPOSPHERE 22 

5. DERIVATIVES OF OBSERVABLES WITH RESPECT TO MODEL PARAMETERS ... 25 

5.1 GEOMETRIC PARAMETERS 25 

5.2 CLOCK PARAMETERS 28 

5.3 TROPOSPHERE PARAMETERS 30 

5.4 SATELLITE PARAMETERS 31 

6. REFERENCES 32 

APPENDICES 

A. Parameter Units and Physical Constants 34 

B. GPSOMC Parameters 35 


TABLES 

I. Periodic Tidally-Induced Variations in UT1 With Periods Less Than 35 Days 8 

II 1980 IAU Theory of Nutation 11 

III. Dependence of the Constants a and b on Tropospheric Model Parameters 23 

B.I. Glossary of GPSOMC Parameters 35 


preceding page blank not filmed 


\\Yi. 





V 


m a n > 


SECTION 1 


INTRODUCTION 

In interpreting measurements of range from satellites of the Global Positioning System (GPS) 
to ground-based receivers, the observables are passed through a multiparameter estimation routine 
(“filter”) to estimate the parameters of a model. This model describes the spacecraft orbits and the 
motions of the Earth-fixed receivers, and supplies to the filter a priori values of computed observables, 
as well as their partial derivatives with respect to model parameters. During 1984-5, software was 
developed at JPL to perform these functions. The purpose of the present report is to describe that 
portion of the software which concerns the modeling of receiver locations, motions of the whole 
Earth, and computation of observables and their partial derivatives. Modeling of satellite orbits 
and parameter estimation form separate units in the software chain and are described in separate 
documents. 

Many aspects of the model required to describe GPS range measurements are identical to the 
modeling developed during the past decade for Very Long Baseline Interferometry (VLB1). Conse- 
quently much of the content of this document, as well as of the associated software package GPSOMC, 
is borrowed from the VLBI modeling and parameter estimation package MASTERFIT and the doc- 
ument giving the details of its models (Sovers and Fanselow, 1987). There are two major differences. 
First, in satellite range measurements, the sources (transmitters) are not at infinity, and they contain 
time standards of their own. Together with the need for orbit determination, this makes the situation 
considerably more complicated than the VLBI case, where signals are received from fixed sources 
at infinite distances. The second, minor, model difference is the lack of necessity of describing the 
passage of the signal through the Earth's ionosphere. It is assumed that ionospheric effects will always 
be removed by performing measurements at two frequencies. 

Three major model components which will be discussed here are: geometry, clocks, and tropo- 
sphere. Section 2 deals with the coordinate frame for the model and establishes methods for calculating 
the position of the receiver in that frame, employing the best current models of whole Earth motions 
and local tidal deformations. This section is nearly identical to the corresponding coordinate frame 
description in the MASTERFIT document, with slight differences reflecting differences in the im- 
plementation of minor effects. The next section (Sec. 3) defines the observables and the intimately 
associated models of behavior of space vehicle and receiver clocks. Section 4 presents the model em- 
ployed to describe the passage of the signal through the troposphere. All the partial derivatives of 
observables with respect to model parameters are given in Sec. 5. Values of the physical constants used 
in the GPSOMC software, as well as a complete list of the parameters presently available for adjustment, 
are given in the appendices. 

As of August 1988, the model description in this report is intended to correspond to the Fortran 
code used to generate the executable GPSOMC.EXE; 114 on the Section 335 LOGOS and XENOS VAX 
computers. This correspondence is believed to hold for the dominant components of the observable 
models. Some of the code implementing seldom-used aspects of modeling (e.g., nutation parameters, 
time-dependent station positions and zenith tropospheres) has not yet been thoroughly checked; it is 
intended that complete testing will be performed in the near future. 
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SECTION 2 


COORDINATE FRAME AND GEOMETRY 


The geometric range is that portion of the distance between a satellite transmitter and an Earth- 
fixed receiver which would be measured by perfect instrumentation, perfectly synchronized, if there 
were a perfect vacuum between the transmitter and receiver. It is by far the largest component of the 
observed range. The main complexity of this portion of the model arises from the numerous coordinate 
transformations necessary to relate the inertial reference frame used for locating the spacecraft to the 
Earth-fixed reference frame in which station locations are represented. 

In the following we will assume that “geocentric inertial reference frame” means a geocentric, 
equatorial frame with the equator and equinox of J2000 as defined by the 1976 IAU conventions 
with the 1980 nutation series (Seidelmann, 1982, and Kaplan, 1981). On the other hand, we will 
mean by “Earth-fixed reference frame” some reference frame tied to the mean surface features of the 
Earth. This is a right-handed version of the CIO reference system with the pole defined by the 1903.0 
pole. Implementation is accomplished by defining the position of one of the fiducial stations, and 
measuring the positions of the other receivers. This section provides the details for the evaluation 
of the geometric range, including receiver coordinates, tidal effects, and the transformation from 
Earth-fixed to geocentric inertial coordinate systems. 


2*1 RECEIVER LOCATION TIME DEPENDENCE 

Normally the receiver position vector r^ 0 and its components in the Earth-fixed reference frame 
r, p , A, z (radius off spin axis, longitude, and height above the equator, respectively) are considered 
to be time-invariant. An alternative formulation introduces the time rates of change of the station 
coordinates as adjustable parameters. The model is linear, with the components of r£ 0 (t) at time t 
expressed as 


r «p = r °. P + r» P (t ~ to) (2.1) 

A = A 0 + A(t - t 0 ) (2.2) 

z = z° + z(t — to) (2.3) 


Here to is a reference epoch, at which the receiver coordinates are (r° p , A 0 , z°). The receiver “position 
vector” may include an antenna phase center offset and an offset to a standard benchmark or monu- 
ment (see Sec. 2.9). In the model development that follows, te 0 includes these offsets, but is referred 
to as the “receiver vector.” Besides the linear time variation of Eqs. (2. 1-2.3), modeling allows for 
estimation of a new benchmark position for every observing session. 


2,2 TIDAL EFFECTS 

In calculating the geometric range, we need to consider the effects of crustal motions on receiver 
locations. Among these deformations are solid Earth tides, tectonic motions, and alterations of the 
Earth’s surface due to local geological, hydrological, and atmospheric processes. In the standard 
Earth-fixed coordinate system, tidal effects modify the receiver location te 0 by an amount 

A = A so i -f* A ocn + A po i (2*4) 

where the three terms are due to solid Earth tides, ocean loading, and pole tide, respectively. Other 
Earth- fixed effects would be incorporated by augmenting the definition of A. 
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2.2.1 Solid Earth Tides 


Alteration of the positions of the receivers by solid Earth tides is rather complicated due to their 
coupling with the ocean tides and to the effects of local geology. We gloss over these complications, and 
employ the simple quadrupole response model described by J. G. Williams (1970), who used Melchior 
(1966) as a reference. Let R* be the position vector of a perturbing source in the Earth- fixed reference 
system, and let Ye 0 be the receiver position vector in the same coordinates. If h and i are the Love 
numbers, rp a phase shift of the tidal effects, and r p the phase-shifted receiver location vector in the 
Earth-fixed reference system, then the vector of tidal displacements in a local Cartesian frame (x axis 
vertical, y axis eastward, and z axis northward on a spherical Earth) is 

* = X]l h 9l*> 192; lgsi] T (2.5) 


where 


9u = 


3 

R* 


(r P R .) 3 _ r l R2 t 

2 6 


3 H.r? 


?2» — £ 5 P ( r p • R»)(^»£p 


K 


*l + y % 


„ _ 3 ^ r p (r R ^ 

93 , ^5~l r P ' J 


yj4 + y 2 P z * - - rf + yp Y >) 


\T x l + yl 


( 2 . 6 ) 

(2.7) 


( 2 . 8 ) 


Here is the ratio of the mass of the disturbing object, s, to the mass of the Earth, and 

R , = [X„Y„Z,} t (2.9) 

is the vector from the center of the Earth to that body. The summation is over all tide-producing 
bodies, of which we include only the Sun and the Moon. 

The phase-shifted receiver vector is calculated employing a phase lag, or equivalently, employing 
a right-handed rotation, L } through an angle rp about the z axis of date, r p = Lr^ 0 . This lag matrix, 
L 1 is: 

( cos rp sin rp 0 A 

— sin^» cos rp 0 I (2.10) 

0 0 ij 

By a positive value of rp we mean that the peak response on an Earth meridian occurs at a time 
St = rp/wE after that meridian containing Te 0 crosses the tide-producing object, where u)e ia the 
angular rotation rate of the Earth. In the vertical component, the peak response occurs when the 
meridian containing r p also includes R,. The difference between geodetic and geocentric latitude can 
affect this model on the order of the (tidal effect)/ (flattening factor) « 0.1 cm. 

To convert the locally referenced strain, S } to the Earth-fixed frame, two additional rotations 
must be performed. The first, W , rotates by an angle, — <p t , about the y axis to an equatorial system. 
The second, V , rotates about the resultant z axis by angle, — A a , to bring the displacements into the 
standard Earth-fixed coordinate system. The result is 


where 


AjoJ = 

= vws 

(2.11) 

/ cos^. 

0 — gin <f> , \ 


3 

ii 

o 

1 ° 

(2.12) 

\ sin <p 0 

0 cos^ J 
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and 


V = 


cos A, 
sin Aj 
0 


(2.13) 



Actually, the product of these two matrices is coded: 

'cos A* cos <f> t — sinA a — cos A, sin <f> 9 

VW — | sin A, cos <p B cos A, — sin A , sin 4> t 

sin <j> 9 0 cos <f> t 


(2.14) 


2.2.2 Ocean Loading 

This section is concerned with one of the secondary tidal effects, t.e., the elastic response of the 
Earth’s crust to ocean tides, which move the receivers to the extent of a few cm. Such effects are 
commonly labeled “ocean loading.” The most complete recent model appears to be that described 
by Pagiatakis (1982, 1988) and by Pagiatakis, Langley, and Vanicek (1982), which is implemented 
as an option in the current GPSOMC. It was recently generalized to include effects of anisotropy and 
viscoelasticity of the Earth; differences of displacements with those given by the 1982 model are well 
below 1 cm. The formulation and coding are general enough, however, to permit other inputs to 
be used in place of the Pagiatakis ocean loading model. Because the receiver motions caused by 
response to ocean tides appear to be limited to approximately 3 cm for sites well removed from the 
coast, no estimation capability was deemed necessary at present. This decision is supported by the 
fact that for locations near the coast, where the effects may be more sizable, and which would thus 
be expected to produce data useful in parameter estimation, the elastic response modeling is as yet 
inadequate (Agnew, 1982). The present model entails deriving an expression for the locally referenced 
displacement 6 due to ocean loading. In the vertical, N-S, E-W local coordinate system at time f, 

N 

Sj = ^ £? cos(w f < + Vi - 6j) (2.15) 

t=i 

The quantities u (frequency of tidal constituent t) and Vi (astronomical argument of constituent t) 
depend only on the ephemeris information (positions of the Sun and Moon). On the other hand the 
amplitude £? and Greenwich phase lag 6/ of each tidal component j are determined by the particular 
model assumed for the deformation of the Earth. As of November 1988, software for calculating 
these deformations at an arbitrary point on the Earth’s surface exists at JPL only for the Pagiatakis- 
Langley model. Six tidal components are included [N — 6 in Eq. (2.15)]: the M 2 , S 2 , -Ki, &U N 2i 
and Pi tides, all of which have periods close to either 12 or 24 hours. The local displacement vector is 
transformed via Eqs. (2.14) and (2.11) to the displacement Aoe» in the standard Earth- fixed frame. 

Input to GPSOMC provides for specification of an arbitrary number of frequencies and astronomical 
arguments and Vi, followed by tables of the local distortions and their phases, and 8* t calculated 
from the ocean tidal loading model of choice. In particular, longer-period tidal constituents can be 
accommodated in this fashion. 

There are presently four choices of models. As mentioned above, the six components of the 
Pagiatakis-Langley model can be calculated by separate software for an arbitrarily located receiver 
to generate input tables for GPSOMC. This has been done for all stations commonly employed in JPL 
VLBI experiments. There is no comparable facility to obtain amplitudes and phases for the Agnew 
(1982), Scherneck (1983), or Goad (1983) models. Consequently for these three models GPSOMC input 
tables only exist for the limited set of stations considered by these authors. Agnew considers only five 
components (omitting P x ), while Scherneck includes five components in addition to the Pagiatakis- 
Langley set: K 2} Q I? M/, M m , and S $a . These all have amplitudes of 1 mm or smaller and, thus, 
are not expected to be significant at the present level of experimental accuracy. Goad’s MERIT 
standard model only specifies vertical displacements, but includes three components (K 2i Q\ and 
Mf) in addition to Pagiatakis-Langley. Due to their bulk, none of the tables of tidal amplitudes are 
reproduced here, but are available on request in computer-readable form. 
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2.2.3 Pole Tide 


In addition to the effects of ocean tides, another secondary tidal effect is the displacement of a 
receiver by the elastic response of the Earth’s crust to shifts in the spin axis orientation. The spin axis 
is known to describe a circle of « 20-m diameter at the north pole. Depending on where the spin axis 
pierces the crust at the instant of a measurement, the “pole tide” displacement will differ from time 
to time. This effect must be included if centimeter accuracy is desired, especially for measurements 
spanning an appreciable fraction of a year. 

Yoder (1984) derived an expression for the displacement of a point at latitude <£, longitude A due 
to the pole tide: 


6 = 


^^[sin<£cos^(a:cos A + y sin A) hr 
9 


+ cos 2 <f>(x cos A + y sin A) 1$ 
+ xsin A -f ycos A)fA] 


(2.16) 


Here a; a is the rotation rate of the Earth, R the radius of the (spherical) Earth, g the acceleration 
due to gravity at the Earth’s surface, and h and l the customary Love numbers. Displacements of the 
spin axis from the 1903.0 CIO pole position along the x and y axes are given by x and y. Eq. (2.16) 

A. 

shows how these map into receiver displacements along the unit vectors in the radial (r), latitude (^), 
and longitude (A) directions. With the standard values ue = 7.292 x 10 -5 rad/sec, R = 6378 km, 
and g = 980.665 g/cm 2 , the factor oj\R/g = 3.459 x 10 “ 3 . Since the maximum values of x and y are 
of the order of 10 meters, and h « 0 . 6 , l « 0.08, the maximum displacement due to the pole tide is 1 
to 2 cm, depending on the location of the receiver ( <f > , A). 

As in the case of the previously considered tidal effects, the locally referenced displacement 6 is 
transformed via the transformation Eq. (2.14) to give the displacement A po t in the standard Earth- 
fixed coordinate system. After each of the locally referenced tidal displacements has been transformed 
to these coordinates, the receiver location is 


** E — ^Eq + A BO i 4* A^cn H- A pol 


(2.17) 


2.3 TRANSFORMATION FROM EARTH-FIXED TO GEOCENTRIC 
INERTIAL COORDINATE SYSTEMS 

The Earth is approximately an oblate spheroid, spinning in the presence of two massive moving 
objects (the Sun and the Moon) which are positioned such that their time-varying gravitational effects 
not only produce tides on the Earth, but also subject it to torques. In addition, the Earth is covered 
by a complicated fluid layer, and also is not perfectly solid internally. As a result, the orientation 
of the Earth is a very complicated function of time, which to first order can be represented as the 
composite of a time-varying rotation rate, a wobble, a nutation, and a precession. The exchange of 
angular momentum between the solid Earth and the fluids on its surface is not readily predictable, 
and thus must be continually determined experimentally. Nutation and precession are well modeled 
theoretically. At the centimeter level, however, even these models are not completely adequate. 

The rotational transformation, Q, of coordinate frames from the Earth-fixed frame to the geocen- 
tric inertial frame is composed of 6 separate rotations (actually 10 , since the nutation transformation, 
N , consists of 3 transformations in itself, as does the “perturbation” transformation, Q) applied to a 
vector in the Earth-fixed system: 

Q = n PNUXY (2.18) 

Thus, if Te is a receiver location expressed in the Earth-fixed system, e.y., the result of Eq. (2.17), 
that location, 17 , expressed in the geocentric inertial (J 2000 ) system is 

17 = Qt e (2.19) 
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Note that since we rotate the Earth rather than the celestial sphere, our rotation matrices, H, P, and 
N y will be the transposes of those used to rotate the inertial system of J2000 to an inertial system of 
date. 

2.4 UT1 AND POLAR MOTION 

The first transformation, Y, is a right-handed rotation about the x axis of the Earth-fixed frame 
by an angle 0 2 . Currently, the Earth- fixed frame is the 1903.0 CIO frame, except that the positive y 
axis is at 90 degrees east (Moscow). The x axis is coincident with the 1903.0 meridian of Greenwich, 
and the z axis is the 1903.0 standard pole. The Y rotation matrix is 

(10 o A 

Y = 0 cos ©2 sin 0 2 I 

\0 —sin © 2 cos © 2 J 

where ©2 is identical to the © y that may be readily obtained from the pole position published by BIH 
(1987) or IRIS (I AG, 1987). 

The next rotation in sequence is the right-hand rotation through an angle ©1 about the y axis 
obtained after the previous rotation has been applied: 

( cos©i 0 — sin©i\ 

0 10 
sin©i 0 cos©i J 

In this rotation, ©1 is identical to the ©* obtainable from the BIH or IRIS value for the pole position. 
Note that we have incorporated in the matrix definitions the transformation from the left-handed 
system used by BIH to the right-handed system we use. Note also that any source of polar motion 
data can be used provided it is represented in a left-handed system. The only effect would be a change 
in the definition of the Earth-fixed reference system. 

The application of u XY n to a vector in the Earth-fixed system of coordinates expresses that 
vector as it would be observed in a coordinate frame whose z axis was along the Earth’s ephemeris 
pole. The third rotation, U , is about the resultant z axis obtained by applying “XY.” It is a rotation 
through the angle, — H, where H is the hour angle of the true equinox of date (t.e., the dihedral angle 
measured westward between the xz plane defined above and the meridian plane containing the true 
equinox of date). The equinox of date is the point defined on the celestial equator by the intersection 
of the mean ecliptic with that equator. It is that intersection where the mean ecliptic rises from below 
the equator to above it (ascending node). 


( 2 . 21 ) 


( 2 . 20 ) 


( cos H —sin H (A 

sin if cos H 0 (2.22) 

001) 

The angle H is composed of two parts: 

H = h 7 + as (2.23) 

where h 7 is the hour angle of the mean equinox of date, and oce is the difference in hour angle of 
the true equinox of date and the mean equinox of date, a difference which is due to the nutation of 
the Earth. This set of definitions is cumbersome and couples the nutation and precession effects into 
Earth rotation measurements. However, in order to provide a direct estimate of conventional UTl it 
is convenient to endure this historical approach, at least for the near future. 

UTl (universal time) is defined to be such that the hour angle of the mean equinox of date is 
given by the following expression (Aoki et al ., 1982, and Kaplan, 1981): 


h n = UTl + 6 h 41 m 50 s . 54841 + 8640184*.812866 T u 

+ 0 s .093104 Tj - 6 s . 2 X 10 -6 T* (2.24) 
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where 


(Julian UT 1 date) - 2451545.0 
u ” 36525 

The actual equivalent expression which is coded is: 

h 7 =2n(UTl Julian day fraction) + 673 10*. 54841 

+ 8640184*. 812866 T u + 0*.093104 T 2 - 6V2 x 10“ 6 T* 


(2.25) 


(2.26) 


This expression produces a time, UT 1, which tracks the Greenwich hour angle of the real Sun to 
within 16™. However, it really is sidereal time, modified to fit our intuitive desire to have the Sun 
directly overhead at noon on the Greenwich meridian. Historically, differences of UTl from a uniform 
measure of time, such as atomic time, have been used in specifying the orientation of the Earth. Note 
that this definition has buried in it the precession constant since it refers to the mean equinox of date. 

By the very definition of “mean of date” and “true of date,” nutation causes a difference in the 
hour angles of the mean equinox of date and the true equinox of date. This difference, called the 
“equation of equinoxes,” is denoted by ole and ls obtained accordingly: 


where the vector 



( X 'A 

* 

U; 



i 

0 

0 


(2.27) 


(2.28) 


is the unit vector, in true equatorial coordinates of date, toward the mean equinox of date. In mean 
equatorial coordinates of date, this same unit vector is just (1,0, 0) T . The matrix is just the 

inverse (or equally, the transpose) of the transformation matrix N which will be defined below to 
effect the transformation from true equatorial coordinates of date to mean equatorial coordinates of 
date. 

Depending on the smoothing used to produce the a priori UTl — UTC series, the short-period 
(t <35 days) fluctuations in UTl due to changes in the latitude and size of the mean tidal bulge may 
or may not be smoothed out. Since we want as accurate an a priori as possible, it may be necessary 
to add this effect to the UTl a priori obtained from the series, UTl 9moot hcd * If this option is selected, 
then the desired a priori UTl is given by 


UTl a priori — UTl 9moo thed “f" AUT1 


(2.29) 


UTl amoothed represents an appropriately smoothed a priori measurement of the orientation of the 
Earth (e.g. } BIH circular D smoothed or, even better, UTlR) y for which the short period (t < 35 days) 
tidal effects either have been averaged to zero, or, as in the case of UT1R, removed before smoothing. 
This AUT1 can be represented as 


N 

■ 

■5 y 

AUTl = ^ 

Ai sin 

X/ k\jOtj 

*=1 


0=1 J. 


(2.30) 


where N is chosen to include all terms with a period less than 35 days. There are no other con- 
tributions until a period of 90 days is reached. However, these long-period terms are included by 
the measurements of the current Earth-orientation measurement services. The values for fc t *y and A*, 
along with the period involved are given in Table I. The for i = 1, 5 are just the angles defined 
below in the nutation series as /, F , D } and 0 respectively. 

It is convenient to apply U UXY ” as a group. To parts in 10 12 , XY = YX. However, with the 
same accuracy UXY ^ XYU . Neglecting terms of O(0 2 ) (which produce receiver location errors of 
approximately 6 x 10“ 6 meters): 


f cos H 
sin H 
\^sin©i 


-sinff 
cos H 
- sin 0 2 


UXY = 


— sin 0i cos H — sin ©2 sin H 

— sin ©1 sin H + sin ©2 cos H 

1 
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(2.31) 


Table I 

Periodic Tidally Induced Variations in UT1 
With Periods Less Than 35 Days 
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2.5 NUTATION 


With the completion of the UT 1 and polar motion transformations, we are left with a receiver 
location vector, r date . This is the receiver location in true equatorial inertial coordinates of date. The 
last set of transformations are nutation, N , precession, P, and the perturbation rotation, O, applied in 
that order. These transformations give the receiver location, rj, in geocentric inertial J2000 equatorial 
coordinates: 

r/ = QPNr da te (2.32) 

Both the nutation and precession rotation angles are defined relative to their values at Julian date 
2451545.0 (J2000). The angles are computed from trigonometric and polynomial series as a function 
of Barycentric Dynamic Time ( TDB ). This time scale, which is also used to reference the planetary 
ephemeris, is related to Terrestrial Dynamic Time ( TDT ) by (Kaplan, 1981) 

TDB = TDT + 0 a . 001658 sin (g + 0.0167sing) (2.33) 

g = 2tt(357°.578 + 35999°.050 TDT) /360° (2.34) 

The time scale TDT runs at the same rate as International Atomic Time (TAI) and is offset from 
TAI by a defined constant, 

TDT = TAI + 32M84 (2.35) 


All time arguments used in nutation and precession computations are measured in centuries of TDB 
from J2000 [cf. Eq. (2.25)]. 

The transformation matrix N is a composite of three separate rotations (Melbourne et a/., 1968): 

1. A(e) : true equatorial coordinates of date to ecliptic coordinates of date. 


A{e) = 



0 

cos e 
— sine 


° 1 

sine 

cos e J 


(2.36) 


where e is the true obliquity of the ecliptic. 

2. C T [Sip) : nutation in longitude from ecliptic coordinates of date to mean ecliptic coordinates of 
date. 


( cos Sip sin Sip 0\ 

— sin Sip cos Sip 0 j (2.37) 

ooi ; 

where Sip is the nutation in ecliptic longitude. 

3. A r (e) : ecliptic coordinates of date to mean equatorial coordinates. 

In ecliptic coordinates of date, the mean equinox is at an angle, Sip = tan -1 ^/^?)* Se = e — e 
is the nutation in obliquity, and e is the mean obliquity (the dihedral angle between the plane of the 
ecliptic and the mean plane of the equator). “Mean” as used in this section implies that the short- 
period (T < 18.6 years) effects of nutation have been removed. Actually, the separation between 
nutation and precession is rather arbitrary, but historical. The composite rotation is: 

N = A T (e)C T {5il>)A{ e ) (2.38) 

( cos Sip cose sin Sip sine sin Sip \ 

— cosesin Sip cosecosecos^t/i + sinesine cos e sin e cos 6 ip — sinecose ) 

— sinesin Sip sinecosecos Sip — cosesine sin esin ecos Sip + cose cos ej 

The 1980 IAU nutation model (Seidelmann, 1982, and Kaplan, 1981) is used for obtaining the 
values for Sip and e — e. The mean obliquity is obtained from Lieske et al. (1977) or from Kaplan 
(1981): 

e = 23° 26' 21."448 - 46 "8150 T - 5."9 X 10 _4 T 2 + l"813 X 10~ 3 T 3 (2.39) 
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This nutation in longitude (6iJ>) and in obliquity (Se = e-e) can be represented by a series expansion 
of the sines and cosines of linear combinations of five fundamental arguments. These are (Kaplan, 
1981, Cannon, 1981): 

1. The mean anomaly of the Moon: 


a x = l= 485866”. 733 + (1325 r + 715922".633) T 

+ 31". 3 10 T 2 + 0".064 T 3 (2.40) 

2. The mean anomaly of the Sun: 

a 2 = l' = 1287099”. 804 + (99 r + 129258l".224) T 

- 0".577 T 2 - 0."012 T 3 (2.41) 

3. The mean argument of latitude of the Moon: 

a 3 = F= 335778". 877 + (1342 r + 295263".137) T 

- 13". 257 T 2 + 0".011 T 3 (2.42) 

4. The mean elongation of the Moon from the Sun: 

a 4 = £> = 107226l".307 + (1236 r + 110560l".328) T 

- 6". 891 T 2 + 0".019 T 3 (2.43) 


5. The mean longitude of the ascending lunar node: 

a 5 = n = 450160".280 - (5 r + 482890".539) T 

+ 7".455 T 2 + 0".008 T 3 (2.44) 


where l r = 360° = 1296000". 

With these fundamental arguments, the nutation quantities then can be represented by 


N 


= £ 


(Aoy + AxjT) sin 


y=il 


.t=i 


(2.45) 


and 


N 


' 6 

* 

Se = J2 

3 = 1 

(Bqj + BijT) cos 

£**<*.• m 

.t=i 



(2.46) 


where the various values of a*, fcy», Ay, and J3y are listed in Table II. 

An additional set of terms can be optionally added to the nutations Srp and 8s in Eqs. (2.45) and 
(2.46). These include the out-of-phase nutations, and the free-core nutations (Yoder, 1983) with period 
ojf (nominally 460 days). The “nutation tweaks” A %p and A s are arbitrary constant increments of 
the nutation angles dtp and 8e . Unlike the usual nutation expressions, they have no time dependence. 
The out-of-phase nutations, which are not included in the IAU 1980 nutation series, are identical to 
Eqs. (2.45) and (2.46), with the replacements sin <-► cos: 


and 


N 6 

^ = ^ > [(-^2y + ^3yT*) c °s ^ ^ A:y»a t (T)j j 

3=1 i= i 

N 5 

Se = E[( 5 2y + 5 3 yT)sin[E^«,(T)]] 
y=i *=i 


(2.47) 


(2.48) 
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Table II 


1 

j 

t 1980 IAU Theory of Nutation 


Index 

j 

Period 

( days ) 

D 

B 

i 


B 

- 4(>y Aij 

( 0 ". 0001 ) 

Bo j Bij 

( 0 ". 0001 ) 

1 

6798.4 

0 

0 

0 

0 

1 

-171996 

- 174.2 

92025 

8.9 

2 

3399.2 

0 

0 

0 

0 

2 

2062 

0.2 

-895 

0.5 

3 

1305.5 

-2 

0 

2 

0 

1 

46 

0.0 

-24 

0.0 

4 

1095.2 

2 

0 

-2 

0 

0 

11 

0.0 

0 

0.0 

5 

1615.7 

-2 

0 

2 

0 

2 

-3 

0.0 

1 

0.0 

6 

3232.9 

1 

-1 

0 

-1 

0 

-3 

0.0 

0 

0.0 

7 

6786.3 

0 

-2 

2 

-2 

1 

-2 

0.0 

1 

0.0 

8 

943.2 

2 

0 

-2 

0 

1 

1 

0.0 

0 

0.0 

9 

182.6 

0 

0 

2 

-2 

2 

-13187 

- 1.6 

5736 

- 3.1 

10 

365.3 

0 

1 

0 

0 

0 

1426 

- 3.4 

54 

- 0.1 

11 

121.7 

0 

1 

2 

-2 


-517 

1.2 

224 

- 0.6 

12 

365.2 

0 

-1 

2 

-2 


217 

* 

o 

bi 

-95 

0.3 

13 

177.8 

0 

0 

2 

-2 


129 

0.1 

1 

o 

0.0 

14 

205.9 

2 

0 

0 

-2 

0 

48 

0.0 

1 

0.0 

15 

173.3 

0 

0 

2 

-2 

0 

-22 

0.0 

0 

0.0 

16 

182.6 

0 

2 

0 

0 

0 

17 

- 0.1 

0 

0.0 

17 

386.0 

0 

1 

0 

0 

1 

-15 

0.0 

9 

0.0 

18 

91.3 

0 

2 

2 

-2 

2 

-16 

0.1 

7 

0.0 

19 

346.6 

0 

-1 

0 

0 

1 

-12 

0.0 

6 

0.0 

20 

199.8 

-2 

0 

0 

2 

1 

-6 

0.0 

3 

0.0 

21 

346.6 

0 

-1 

2 

-2 

1 

-5 

0.0 

3 

0.0 

22 

212.3 

2 

0 

0 


1 

4 

0.0 

-2 

0.0 

23 

119.6 

0 

1 

2 

-2 

i 

4 

0.0 

-2 

0.0 

24 

411.8 

1 

0 

0 

-1 


-4 

0.0 

0 

0.0 

25 

131.7 

2 

1 

0 

-2 

0 

1 

0.0 

0 

0.0 

26 

169.0 

0 

0 

-2 

2 

1 

1 

0.0 

0 

0.0 

27 

329.8 

0 

l 

-2 

2 

0 

-1 

0.0 

0 

0.0 

28 

409.2 

0 

l 

0 

0 

2 

1 

0.0 

0 

0.0 

29 

388.3 

-1 

0 

0 

1 

1 

1 

0.0 

0 

0.0 

30 

117.5 

0 

1 

2 

-2 

0 

-1 

0.0 

0 

0.0 

31 

13.7 

0 

0 

2 

0 

2 

-2274 

- 0.2 

977 

- 0.5 

32 

27.6 

1 

0 

0 

0 

0 

712 

0.1 

-7 

0.0 

33 

13.6 

0 

0 

2 

0 

1 

-386 

- 0.4 

200 

0.0 

34 

9.1 

1 

0 

2 

0 

2 

-301 

0.0 

129 

- 0.1 

35 

31.8 

1 

0 

0 

-2 

0 

-158 

0.0 

-1 

0.0 

36 

27.1 

-1 

0 

2 

0 

2 

123 

0.0 

-53 

0.0 

37 

14.8 

0 

0 

0 

2 

0 

63 

0.0 

-2 

0.0 

38 

27.7 

1 

0 

0 

0 

1 

63 

0.1 

-33 

0.0 

39 

27.4 

-1 

0 

0 

0 

1 

-58 

- 0.1 

32 

0.0 

40 

9.6 

-1 

0 

2 

2 

2 

-59 

0.0 

26 

0.0 
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Table II cont. 


1980 IAU Theory of Nutation 


Index 

j 

Period 

(days) 

kji 


kj3 


kjs 

•doy ^iy 

(0".0001) 

Boj Bij 

(0".000l) 

41 

9.1 

1 

0 

2 

0 

1 

-51 

0.0 

27 

0.0 

42 

7.1 

0 

0 

2 

2 

2 

-38 

0.0 

16 

0.0 

43 

13.8 

2 

0 

0 

0 

0 

29 

0.0 

-1 

0.0 

44 

23.9 

1 

0 

2 

-2 

2 

29 

0.0 

-12 

0.0 

45 

6.9 

2 

0 

2 

0 

2 

-31 

0.0 

13 

0.0 

46 

13.6 

0 

0 

2 

0 

0 

26 

0.0 

-1 

0.0 

47 

27.0 

-1 

0 

2 

0 

1 

21 

0.0 

-10 

0.0 

48 

32.0 

-1 

0 

0 

2 

1 

16 

0.0 

-8 

0.0 

49 

31.7 

1 

0 

0 

-2 

1 

-13 

0.0 

7 

0.0 

50 

9.5 

-1 

0 

2 

2 

1 

-10 

0.0 

5 

0.0 

51 

34.8 

1 

1 

0 

-2 

0 

-7 

0.0 

0 

0.0 

52 

13.2 

0 

1 

2 

0 

2 

7 

0.0 

-3 

0.0 

53 

14.2 

0 

-1 

2 

0 

2 

-7 

0.0 

3 

0.0 

54 

5.6 

1 

0 

2 

2 

2 

•8 

0.0 

3 

0.0 

55 

9.6 

1 

0 

0 

2 

0 

6 

0.0 

0 

0.0 

56 

12.8 

2 

0 

2 

-2 

2 

6 

0.0 

-3 

0.0 

57 

14.8 

0 

0 

0 

2 

1 

-6 

0.0 

3 

0.0 

58 

7.1 

0 

0 

2 

2 

1 

-7 

0.0 

3 

0.0 

59 

23.9 

1 

0 

2 

-2 

1 

6 

0.0 

-3 

0.0 

60 

14.7 

0 

0 

0 

-2 

1 

-5 

0.0 

3 

0.0 

61 

29.8 

1 

-1 

0 

0 

0 

5 

0.0 

0 

0.0 

62 

6.9 

2 

0 

2 

0 

1 

-5 

0.0 

3 

0.0 

63 

15.4 

0 

1 

0 

-2 

0 

-4 

0.0 

0 

0.0 

64 

26.9 

1 

0 

-2 

0 

0 

4 

0.0 

0 

0.0 

65 

29.5 

0 

0 

0 

1 

0 

-4 

0.0 

0 

0.0 

66 

25.6 

1 

1 

0 

0 

0 

-3 

0.0 

0 

0.0 

67 

9.1 

1 

0 

2 

0 

0 

3 

0.0 

0 

0.0 

68 

9.4 

1 

-1 

2 

0 

2 

-3 

0.0 

1 

0.0 

69 

9.8 

-1 

-1 

2 

2 

2 

-3 

0.0 

1 

0.0 

70 

13.7 

-2 

0 

0 

0 

1 

-2 

0.0 

1 

0.0 

71 

5.5 

3 

0 

2 

0 

2 

-3 

0.0 

1 

0.0 

72 

7.2 

0 

-1 

2 

2 

2 

-3 

0.0 

1 

0.0 

73 

8.9 

1 

1 

2 

0 

2 

2 

0.0 

-1 

0.0 

74 

32.6 

-1 

0 

2 

-2 

1 

-2 

0.0 

1 

0.0 

75 

13.8 

2 

0 

0 

0 

1 

2 

0.0 

-1 

0.0 

76 

27.8 

1 

0 

0 

0 

2 

-2 

0.0 

1 

0.0 

77 

9.2 

3 

0 

0 

0 

0 

2 

0.0 

0 

0.0 

78 

9.3 

0 

0 

2 

1 

2 

2 

0.0 

-1 

0.0 

79 

27.3 

-1 

0 

0 

0 

2 

1 

0.0 

-1 

0.0 

80 

10.1 

1 

0 

0 

-4 

0 

-1 

0.0 

0 

0.0 
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It is emphasized that, for the present, the default nutation model in GPSOMC is just the 1980 IAU 
nutation model. 


2.6 PRECESSION 

The next transformation in going from the Earth-fixed frame to the geocentric inertial frame is 
the rotation P. This is the precession transformation from mean equatorial coordinates of date to the 
equatorial coordinates of the reference epoch ( e.g.> J2000). It is the transpose of the matrix given on 
page 7 of Melbourne et al. (1968): 


P lt = cos cos ©x cos Za — sin £4 sin Za (2.53) 

P 12 = cos £4 cos sin Za 4- sin £4 cos Za (2.54) 

P 13 = cos ^ sin 0 a (2.55) 

P 21 = — sin $ a cos &a cos Za — cos £4 sin Za (2.56) 

P 22 = — sin $a cos ©x sin Za 4- cos £4 cos Za (2.57) 

P 23 = — sinfxsin©x (2.58) 

P 31 = — sin ©x cos Z A (2.59) 

P 32 = — sin ©x sin Z A (2.60) 

P 33 = cos ©>i (2.61) 

With the angular units in arc seconds, the arguments are: 

fx = 2306.2181 T + 0.30188 T 2 + 0.017998 T 3 (2.62) 

Z A = 2306.2181 T + 1.09468 T 2 + 0.018203 T 3 (2.63) 

©x = 2004.3109 T - 0.42665 T 2 - 0.041833 T 3 (2.64) 


These expressions are given by Lieske et al. (1977) and by Kaplan (1981). This completes the standard 
model for the orientation of the Barth. 
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2.7 PERTURBATION ROTATION 


The standard model for the rotation of the Earth as a whole may need a small incremental 
rotation about any one of the resulting axes, for example, in compensating for defects in the a priori 
precession model. Define this perturbation rotation matrix as 

A X A„A, (2.65) 


0 0 \ 

1 -50, ( 2 . 66 ) 

se x 1 J 

with 50, being a small angle rotation about the x axis, in the sense of carrying y into z; 

( i o se„\ 

A„ = 0 10 (2.67) 

V-50„ 0 1 J 

with 5© y being a small angle rotation about the y axis, in the sense of carrying z into x; and 

( 1 -50, 0 ^ 

50, 1 0 ( 2 . 68 ) 

0 0 l) 

with 50, being a small angle rotation about the z axis, in the sense of carrying x into y. For angles 
of the order of 1 arc second we can neglect terms of order SB 2 Re as; they give effects on the order 
0.015 cm. Thus, in that approximation 


where 


n 


A* = 


f 1 - 50 , 6 By \ 

O = 50, 1 -6B X (2.69) 

V -5e y SB X 1 J 

In general, 

50; = 50<(t) = 50 io +50 i T+/ < (r) (2.70) 

which is the sum of an offset, a time-linear rate, and some higher-order or oscillatory terms. Currently, 
only the offset and linear rate are implemented. In particular, 50„ is equivalent to a change in the 
precession constant. Setting 

50, = 6B y = 50, = 0 (2.71) 

gives the effect of applying only the standard rotation matrices. 

Starting with the Earth-fixed vector, Te 01 we have in Sections 2 . 1 — 2.6 above shown how we obtain 
the same vector, 17 , expressed in the geocentric inertial frame: 

rj = nPNUXY{r Eo + A) (2.72) 
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2.8 GEOCENTER OFFSET AND COORDINATE SCALING 


The Earth-fixed reference frame is essentially derived from VLBI measurements, which are com- 
pletely insensitive to the location of the Earth’s center of mass. Practically, receiver coordinates are 
based on the location of a reference station, which is derived either from spacecraft tracking data or 
satellite laser ranging (SLR) measurements. Consequently, there might be sizable errors in fiducial 
station locations (~ 10 m or ~ 1 m for the two techniques, respectively). GPS tracking data is capable 
of determining the position of the geocenter to an accuracy dependent on the quality and quantity of 
the range observations. In order to allow for the possibility of solving for the location of the geocenter, 
Cartesian offsets were introduced, such that the receiver coordinates (including those of the reference 
station(s)) are 


X E 0 X E 0 + X GC 

VEq — ► Ve 0 + Vgc (2.73) 

ze 0 —* ZEo + zqc 

or in vector notation, 

TEo r Eo + r GC (2.74) 

Another problem with coordinate systems at the present level of understanding is the uncer- 
tainty that the coordinate transformations for a priori satellite orbits are being carried out correctly. 
Introduction of a scale factor a, which multiplies Earth-fixed coordinates relative to the geocentric 
inertial coordinates of the satellites, permits empirical detection of such problems. It gives rise to the 
transformations 


*E 0 ax E 0 


yEo -*• <*ys 0 

(2.75) 

z Eo - *■ az E 0 


TBo - *■ ar Eo 

(2.76) 


in vector notation, for the coordinates of the receivers. 

Upon inclusion of both of these modifications in the Earth-fixed station locations Te 0 of Eq. 
(2.72), the final expression for the station location vector in geocentric inertial coordinates becomes 

rj = ClPNUXY(a r Eo 4- r G c + A) = Q[oite 0 4- r G c + A) (2.77) 


2.9 PHASE CENTER OFFSETS 

Geometric offsets may exist between transmitter and receiver phase centers and the standard 
reference points. In the case of the ground-based receiver, this offset takes the form of a “site vector” 
from a surveyor’s benchmark to the receiver phase center. Normally, the position of the benchmark, 
*bm, is desired, rather than the position of the phase center, Tpc . These two vectors are related 
through the site vector Tsurv : 


te 0 = t bm = Tpc 4- tsurv (2.78) 

For GPS satellites, the ephemeris reference point is the spacecraft center of mass (CM). The vector 
offset between the phase center and CM is given by 0.211 i 4-0.886 k meters in a coordinate system in 
which the unit vector k points from the spacecraft CM to the Earth center, j is the normalized cross 
product of k with the unit vector from the spacecraft to the Sun, and i completes a right-handed 
system (Winn, 1984). Both of these geometric offsets are incorporated into GPSOMC modeling. 
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SECTION 3 


OBSERVABLES AND CLOCK PARAMETERS 


Range measurements from GPS satellites are based on the detection of the phase of transmitted 
electromagnetic signals. The physical observable is the difference between the reference phase of the 
transmitter at the time of signal emission and the reference phase of the receiver at the time of signal 
reception. Such measurements of necessity involve detailed and careful consideration of a number of 
time scales. These are intimately related to the definitions of the observables, and for that reason this 
section contains detailed descriptions of both the observables and clock models. 

The GPS satellite transmissions are one-way transmissions originating at the spacecraft. A carrier 
signal modulated by a pseudo-random noise code is broadcast (Spilker, 1978). Both a “carrier phase” 
observable <p and a “pseudo-range” observable £ are modeled inGPSGMC. The observed observables are 
related in a prescribed way to the physically detected phase. The computed observables are defined 
in terms of clock differences since the reference phase at the receiver is derived from the station clock 
and the reference phase at the transmitter is derived from the space vehicle clock. For both carrier 
phase and pseudo-range the computed observable is given by the theoretical difference between the 
space vehicle and station clocks plus a bias term. The clock difference is accumulated as the sum of 
five components — the first is the difference between station clock time and proper time at the station; 
the second is the difference between proper time and coordinate time at the station; the third is the 
difference between the coordinate time of signal emission and the coordinate time of signal reception; 
the fourth is the difference between coordinate time and proper time at the spacecraft; the fifth is the 
difference between proper time at the spacecraft and spacecraft clock time. 

These terms are readily interpreted. The first term is referred to as a station clock error and 
the fifth term is referred to as a space vehicle clock error. The second and fourth terms are general 
relativistic time transformations. The third term is the solution to the light time equation connecting 
the events of signal emission and signal reception. Of these five contributions the third term is 
normally dominant. This term is also referred to as the geometric range between spacecraft and 
receiver, whence the use of “range” in referring to GPS observables. Tropospheric and ionospheric 
delays are not considered here; they are discussed in Sections 1 and 4 of this report. The observables 
are assumed to have been calibrated to remove the effects of instrumental delays. A detailed model 
for instrumental phase delay is not provided, but a bias term is included in the definition of both data 
types. This bias might account for 

(i) an unknown integer number of cycles of phase, 

(ii) an uncalibrated offset in the absolute phase of an oscillator, 

(iii) uncalibrated delays within transmitter or receiver electronics, or 

(iv) an uncalibrated offset between the phase reference used for signal detection and the phase 
reference used for time-tag generation within a receiver. 

The geometric portion of the observable is calculated in the geocentric reference frame, moving 
with the Earth, with inertial J 2000 oriented axes. This reference frame is the generalization of the 
local inertial reference frame containing the Earth (Ashby and Bertotti, 1984) and provides a simple 
setting for applying all necessary relativistic corrections. Terrestrial Dynamic Time (TDT) is used as 
coordinate time (Thomas, 1975 and Ashby and Allan, 1979). The spacecraft equations of motion are 
integrated in this reference frame using TDT as the independent variable. As described in Section 
2 of this report, station coordinates are transformed into this reference system. Care must be taken 
when selecting a priori values for model parameters since particle mass and physical distance in 
the geocentric system differ by a scale factor from current determinations of same in the celestial 
coordinate system (Hellings, 1986). 

We use *2 denote the time of signal emission and £ 3 to denote the time of signal reception. 
Let £ 3 , f 3 , and t 3 refer to, respectively, station clock time, proper time, and coordinate time at the 
receiver. Let £ 2 ) ? 2 > and £3 refer to, respectively, space vehicle clock time, proper time, and coordinate 
time at the transmitter. 
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3.1 PSEUDO-RANGE 


Observed pseudo-range is given by 


Z° = c(h-t 2 ) (3.1) 

with time tag £ 3 , where the times £3 and t 2 are the actual times kept by the station and spacecraft 
clocks, respectively. Correlation of the received pseudo-random noise code generated and transmitted 
by the space vehicle with a local copy of the code within the receiver allows for direct conversion 
between the detected received phase and clock readings. 

Computed pseudo-range is given by 

R c = c ( £3 — £2 + Bstn ) (3.2) 

with time tag £ 3) where now the times £3 and t 2 are the modeled times of, respectively, the station 
and space vehicle clocks. To facilitate computation, this expression is rewritten as a sum of six terms, 

Z C = C ( £3 - £3 
+ £3 ~ fs 
+ £3 “ ^2 
+ t 2 — £2 
+ £2 “ £2 

+ Bstn ) (3.3) 

clock differences discussed above, and the last term is the bias discussed 
modeled as a quadratic function, 

aSTN.Z + bsTN,Jl{t3 - t$ 0 ) + c$tn,jz (^3 ^3 0 ) (3.4) 

a specified epoch. The coefficients asTN y z 1 bsTN,z, and cstn.Z are 
generally unknown but presumed to be small. Ideally, in geocentric coordinates, all clocks fixed to 
the surface of the Earth run at approximately the same rate, with the precise rate depending weakly 
on the geodetic coordinates of the clock. We assume that all clocks are adjusted as necessary so that 
their average rate is the same as that of coordinate time. The rate of coordinate time has been defined 
to agree with the SI second (Kaplan, 1981). Thus, in our model, station proper time, ideal UTC time, 
and coordinate time all run at the same rate. This is not a limitation since the deviation of the rate 
of any clock from the assumed rate can be modeled as a station clock error. 

The difference between proper time and coordinate time at the station is given by 

£3 - * 3 = - [ {TAI - UTC) + (TDT - TAI) ] (3.5) 

where TAI — UTC is an integer number of leap seconds which changes approximately once a year 
and TDT — TAI is defined to be 32.184 sec. 

The geometric range is given by the solution to the light time equation, 

= + Ai„, (3.6) 

c 

where rsrjv^s) is the inertial receiver position at the time of signal reception, rs^ (t 2 ) i s the space 
vehicle position at the time of signal transmission, and A t re i geom is the general relativistic correction to 
the geometric range. To begin the computation of the range, we start from the measurement epoch £ 3 . 
This is converted to £3 using Eq. (3.4) and to £3 using Eq. (3.5). The station position in geocentric 
inertial coordinates is then calculated at the epoch £3 using Eq. (2.77). Thus station clock errors 


The first five terms are the five 
above. 

The station clock error is 
£3 — £3 = 

On the right hand side £ 3o is 


18 



affect the computation of the geometric range. Given the station position and a PV file containing 
the geocentric inertial spacecraft position as a function of coordinate time, the light time equation is 
iteratively solved to obtain t 2 \ 


t 0+i) 

Z 2 




*3 “ 4 f) r 12 / c - 


r « 

"l2 


> W /("S) 


(3.7) 


The space vehicle position used in this computation is that of the antenna phase center, as discussed in 
Section 2.9 of this report. The gravitational delay At re i )eom is given by (Tausner, 1966 and Holdridge, 
1967) 


_ (1 +' 1 ppn)he ,_ri + r2 + ri2 

relgeom 3 . 

c° ri -+• r 2 — ri 2 


(3.8) 


where 7 ppjv is 1 for general relativity, /ip is the Earth’s gravitational constant, and, as in Eq. (3.7), 


r 2 = r sv (* 2 ) 
ri 2 = ri - r 2 


ri = In I 
n = |n I 
r\2 = |ri2| 


(3.9) 


We use rs^fe) as the initial estimate for Tsv{t 2 )\ the light time iteration converges very fast. 

The rate of an ideal space vehicle clock will differ from the rate of coordinate time due to the 
difference in gravitational potential and speed between the space vehicle clock and the reference clock. 
This rate difference may be decomposed into two components: a bias component which depends only 
on the nominal semi-major orbit axis of the space vehicle, and a periodic component which depends 
mainly on the orbit eccentricity. The frequencies of the GPS satellite clocks were purposely set 
low before launch, relative to the nominal published values, to offset the nominal bias rate difference 
between space vehicle time and coordinate time (Spilker, 1978). Accordingly, our model rate difference 
between proper time and coordinate time at the spacecraft contains only a periodic component. This 
component depends only on orbit eccentricity. A bias component due to off-nominal semi-major orbit 
axis and a small periodic component due to Earth oblateness and solar perturbations are absorbed in 
the space vehicle clock error model. The difference between proper time and coordinate time at the 
space vehicle is given by 


t 2 - h = {TAI - UTC) + {TDT - TAI) + 


(3.10) 


where the clock correction term is 


A trel etk = 2 rsv(t 2 ) *rsv(* 2 ) 


(3.11) 


This result may be derived from integration of the usual differential equation relating coordinate 
and proper time [e.y., Eq. (1) of Thomas (1975)] in the geocentric reference system, modeling the 
space vehicle orbit about a point-mass Earth as elliptical while ignoring the differential gravitational 
potential due to perturbing bodies. The resulting rate expression is accurate to about one part in 10 12 
for the GPS orbits, which is comparable to the instability of the GPS satellite clocks. Note that the 
large (« 50 sec) defined offset between coordinate and proper time appears in Eqs. (3.5) and (3.10) 
with opposite sign, and hence does not affect the computation. It is necessary, though, to include this 
offset due to our convention of basing observable time tags on UTC while indexing the PV file by 
TDT. 

Analogously to the station clock error, the space vehicle clock error is also modeled as a quadratic 
function, 

*2 — *2 = a$v,R + bsv,z(t2 — ^2o) + csv,x (*2 “ t 2o ) 2 (3.12) 
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On the right hand side £ 2o is a specified epoch. The space vehicle clock error does not affect geometric 
range. The last term in the model is the pseudo-range bias BstNj which depends on the station but 
is independent of the spacecraft. It also does not affect the value of the computed geometric range. 

3.2 CARRIER PHASE 

Observed carrier phase is given by 

<P° = (-c/w n ) {<f>sv “ <Pstn) (3.13) 

with time tag £3, where <f>sv is the received RF (radio frequency) phase of the space vehicle signal 
at time £3, < f>sTN is the receiver reference phase at time £3, and u) n is the nominal value for both the 
transmitted frequency of the space vehicle signal and the receiver mixing frequency. 

The receiver reference phase is modeled as 


(£3) = W ft ( £3 ” *3 0 ) 


(3.14) 


Since phase is a physical invariant, the received phase of the space vehicle signal at the time of signal 
reception is the same as the space vehicle reference phase at the time of signal transmission, 

<f>Sv{t3) = <t>ref{t2) (3.15) 

The space vehicle reference phase is modeled as 

2 ) = ( $ 2 ” ( 3 . 16 ) 

Thus the computed carrier phase is given by 

<p° = c + Bsv.stn) ( 3 . 17 ) 

with time tag £3. The integration constants t' 3o and t f 2o have been absorbed in the carrier phase bias 
BsVjSTN • Deviation of the station mixing frequency from nominal is modeled as a station clock error, 
and deviation of the space vehicle transmitter frequency from nominal is modeled as a space vehicle 
clock error. Note that the modeled value for carrier phase appears to be the same as the modeled 
value for pseudo-range, except for the substitution of the carrier phase bias for the pseudo-range bias. 
A further distinction is that we allow the modeled values of station time and spacecraft time to be 
different for the two data types. 

Just as for pseudo-range the expression for computed carrier phase is rewritten as a sum of six 
terms, 


<p c = 


c ( £3 - £3 
“f £3 — *3 
+ £3 ~ £2 
+ £2 £2 

+ £2 “ £2 
+ BsVtSTN ) 


The station clock error is modeled as 

£3 “ £3 = aSTN y <p + bsTN,<p{l 3 “ £ 3 o) + c 5 TiV,^ (£3 ~ £30 ) 2 


(3.18) 


(3.19) 


20 



and the space vehicle clock error is modeled as 


*2 — *2 = a SV,v + bsv t <p{t2 — *2 0 ) + C SV,<P (*2 “ *2 0 ) 2 (3.20) 

The carrier phase bias Bsv,stn depends on both the station and the space vehicle. The other terms 
in Eq. (3.18) are evaluated just as for pseudo-range. They will, however, have different numerical 
values unless the station clock error coefficients are identical for the two data types. The value of the 
station clock error affects the geometric range, while the values of the carrier phase bias and space 
vehicle clock error do not. 

In cases where the behavior of transmitter and receiver clocks permits, both the range and phase 
observables may be modeled with a common set of parameters. By analogy with Eqs. (3.4), (3.12), 
(3.19), and (3.20), the common clock errors may be written as 

*3 — *3 = &STN + ^STJV(^3 “ fs 0 ) + C STN (*S “ *3 0 ) (3.21) 

and 

f 2 — f 2 = asv + bsvih — tio) + c$v {fa — t 2 0 ) 2 (3.22) 
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SECTION 4 


TROPOSPHERE 

The delay experienced by the incoming signal due to the Earth’s atmosphere can be modeled 
using a spherical-shell troposphere consisting of a wet component and a dry component. If E is the 
apparent geodetic elevation angle of the spacecraft, the tropospheric contribution to the range is 

Ptrop = PZ ir yRdry{E) + P Z wet Rwct(E) (4.1) 

where pz is the additional delay at local zenith due to the presence of the troposphere, and R is an 
elevation angle mapping function. For recent VLBI data, it was found that modeling the zenith delay 
as a linear function of time improves troposphere modeling considerably. The dry and wet zenith 
parameters are therefore written as 


PZi, w = P°Z iiW + fai.v (t - <o) 


(4.2) 


where to is a reference time. 

The analytic mapping function developed by Lanyi (1984) is used for mapping zenith values to 
the observed elevation angles. In its simplest form, this mapping function employs average values of 
atmospheric constants. Provision is made for specifying surface meteorological data acquired at the 
time of the experiments, which may override the average values. An approximate partial derivative is 
obtained with respect to one parameter in the Lanyi mapping function; this permits adjustment even 
in the absence of surface data. 

Here we attempt to give a minimal summary of the final formulas. The tropospheric delay is 
written as: 

Ptrop = F(E )/ sin E (4.3) 

where 


F(E) = PZ irit Fdry{E) + PZ wei Fwet{E) 

+ [p 2 Zdry F b i(E) + 2 P z dry Pz„ME) + pl wet F b 3 (E)]/A + p% dry F u (E)/A 2 (4.4) 

The quantities pz iry an ^ Pz W9t have the usual meaning: zenith dry and wet tropospheric delays. A 
is the atmospheric scale height, A = kTo/mg C} with k = Boltzmann’s constant, To = average surface 
temperature, m = mean molecular mass of dry air, and g c ~ gravitational acceleration at the center of 
gravity of the air column. With the standard values k = 1.38066 x 10“ 16 erg/K, m = 4.8097 x 10“ 23 
g, g c = 978.37 erg/g cm, and the average temperature for DSN stations To = 292 K, the scale height 
A = 8567 m. 

The dry, wet, and bending contributions to the delay, Fd ry (E) i F wet (E), and F bx , 62 , 63 , 64 (#), axe 
expressed in terms of moments of the refractivity as 


Edry{E) = Ai 0 (i£)G(AAfiio, u) + 3auM2ioG 3 (Miio,u)/4 (4.5) 

F we t{E) = Aoi{E)G(\Mioi/Mooi) u)/Afooi (4.6) 

F b i(E) — [<tG 3 (Mi\o 1 tt)/sin 2 E — Mo 2 o^ 3 (A/i 2 o/Afo 20 j u )]/2 tan 2 E (4*7) 

F b 2{E) = —MoiiG 3 (Mm/Moii ) u)/2Mooitdin 2 E (4*8) 

F b 3 (E) = -Moo2G 3 (M 10 2/Moo2, u)/2 M 0 2 01 tan 2 E (4.9) 

F b t{E) ~ Mo 3 oG 3 (Miso/Mq 3 0 ) u )/2 tan 4 E (4*10) 


A misprinted sign in the last of Eqs. (5) of Appendix B of Lanyi (1984) has been corrected in Eq. 
(4.10). Here G{q>u) is a geometric factor given by 

G{q,v) = (l + 9 u) -1 / 2 (4.11) 
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with 


u = 2a/ tan 2 E (4.12) 

where cr is a measure of the curvature of the Earth’s surface A /R with standard value 0.001345. 

The quantities Ai^E) and are related to moments of the atmospheric refractivity, and are 

defined below. Axo(E) involves the dry refractivity, while A 0 i(2?) is the corresponding wet quantity. 
The Aim{E) are given by 


10 n 

Aim{E) = Molm + EE 


n=lfc=0 


(— l) n+fc (2n — l)!!A/„_fc,i, m 

u 

n 

AA/ifm 

2"fc!(n-Jfc)! 

1 + XuM Um /Moim 


Molm 


(4.13) 


with the scale factor A = 3 for E < 10° and A= 1 for E > 10°. Only the two combinations {l, m) = 
(0,1) and (1,0) are needed for the At m (E). The moments of the dry and wet refractivities are defined 
as 

OO 

M ni} = J q n f' d ryMfLMdq (4.14) 

0 

where fdry t wet{q) are the surface-normalized refractivities. Here, n ranges from 0 to 1, i from 0 to 3, 
and j from 0 to 2; not all combinations are needed. Carrying out the integration in Eq. (4.14) for a 
three-section temperature profile gives an expression for the general moment M n ij : 

AWn! = (1 - c-“^)/« n+1 + [l - n +n+1 (quq?.)} ft 

i=0 

+ e~ aqi T% +n+1 (qi, q 2 )/ a n+1 (4.15) 


Here, 


T 2 (qi,q2) = l- (q2-qi)/<* 


(4.16) 


The quantities q\ and q 2 are the scale-height normalized inversion and tropopause altitudes, respec- 
tively. For the standard atmospheric model, gi = 0.1459 and q 2 = 1.424. The constants a and b are 
functions of the dry (a = 5.0) and wet (/? = 3.5) model parameters, as well as of the powers of the 
refractivities (t and y) in the moment definitions. Table III gives the necessary a’s and 6*s. 


Table III 

Dependence of the Constants a and 6 
on Tropospheric Model Parameters 


B 

m 

a 

b 

1 

0 

1 

a-1 

0 

1 


a/3 — 2 

2 

0 

2 

2(a-l) 

1 

1 

/9+1 

p(a + 1) - 3 

0 

2 

2^ 

2 (ctP - 2) 

3 

0 

3 

3(a-l) 


Note that the normalization is such that Afoio = 1; this moment has therefore not been explicitly 
written in Eqs. (4.5)— (4.10). 
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At present, provision is made for input of four meteorological parameters to override the standard 
(average) values of the Lanyi model. These are: (l) the surface temperature To, which determines the 
atmosphere scale height (default value 292 K); (2) the temperature lapse rate W , which determines 
the dry model parameter a (default values W — 6.8165 K/km, a = 5.0); (3) the inversion altitude 
hi t which determines qi — hi/A (default value hi — 1.25 km); and (4) the tropopause altitude / 12 , 
which determines q 2 = / 12 /A (default value /12 = 12.2 km). A fifth parameter, the surface pressure 
p 0 , is not used at present. Approximate sensitivity of the tropospheric delay (at 15° elevation) to the 
meteorological parameters is —0.3 cm/K for surface temperature, 1 cm/(K/km) for lapse rate, and 
— 1 cm/km for inversion and 0.2 cm/km for tropopause altitude, respectively. 


24 


SECTION 5 


DERIVATIVES OF OBSERVABLES WITH RESPECT TO MODEL PARAMETERS 


This section presents expressions for the partial derivatives of the computed observables with 
respect to the three classes of model parameters described in this document, as well as with respect 
to dynamic parameters related to the spacecraft state. The relativistic corrections, the tropospheric 
correction, and the space vehicle clock error term in the computed observable depend weakly on 
geometric parameters, on the station clock error, and on dynamic parameters. These functional 
dependences are neglected when computing partial derivatives. 

We begin with some notation. From Eqs. (3.2) through (3.17) the computed observable can be 
written as 

O — c ( £3 — £ 3 ) + p + c ( 1 2 — h) + c B u (5*1) 

where C is HP for pseudo-range and yP for carrier phase, and 1 / is STN for pseudo-range and 
5V, STN for carrier phase. 

The clock error terms are given by 

*3“ *3 = &STN,p + t>STN,p(t 3“ *3 0 ) + C STN,P (*3 *3 0 ) 

and 

i 2 — t 2 = asv,p + bsv,p(t 2 — ho) + c sv,p(h — t 2o ) 2 ( 5 - 3 ) 

where /? is H for pseudo-range and <p for carrier phase. The geometric range is given by 

p = | r stn — *sv I ( 5 . 4 ) 

where tstn is the geocentric inertial station position at the time *3 of signal reception and r$v is 
the space vehicle position, in the same system, at the time 1 2 of signal emission. The station position 
comes from Eq. (2.77), 

tstn — fiPNUXY ( arsTNo + t gc + A ) ( 5 . 5 ) 

The space vehicle position and the times of signal emission and reception are related through the 
light time equation, from which we may write the geometric range as 

p = c ( t 3 - t 2 ) 

and, equivalently, 

c 2 (*3 t 2 ) 2 = ( tstn - * sv ) * ( tstn - t S v) ( 5 . 7 ) 



( 5 . 2 ) 


5.1 GEOMETRIC PARAMETERS 


The class of “geometric* parameters includes receiver coordinates and their rates, Love numbers 
and tide phase, relativistic 7 , geocenter offsets, coordinate scale factor, UTPM parameters, nutation 
amplitudes, and the coefficients of the perturbation rotation. Computed observables for both the 
carrier phase and pseudo-range have identical functional dependence on these geometric parameters 
through the geometric range. The partials will therefore be written for the geometric range given by 
Eq. ( 5 . 4 ). Symbolizing one of the geometric parameters by tj, we have 


dp_ 

dr\ 


— -(tstn — * sv) T -£-(tsTN 
p dr\ 




1 . \ dr stn drsv dr stn] 

— ( tstn r s^) ~ 


= -(tstn — Tsv) T 
p 


[*- 


drj 9 tstn 

drsv 1 9tstn 


dTSTN J dl] 


( 5 . 8 ) 
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The light time correction term is computed first: 


drsv . 2 

’ — Y SV 

3rsr/y drsrjv 


(5.9) 


The partial on the right hand side is obtained from Eq . (5.7) by implicit differentiation: 
2c 2 (t 3 - t 2 ) (~g^) = 2(rs™ - r S v) T [/» - 


(5.10) 


Solving for 


dt 2 

drsTN 


gives 


dt 2 

d*STN 


— (ystn — 


cp 




*STN ~ *SV 


r -f) 


(5.11) 


The computation of the geometric partials is complete except for the parameter-dependent terms 
drsTN 
dri 


Taking the parameters in order of their appearance on the right hand side of Eq. (5.5), the partial 
derivatives are 

OYSTN _ dfl -PNUXY(aYsTNo + t gc + A ) ( 5 . 12 ) 


a se 




d se 


*. v,*o 


drsTN 

d 


an 


a se.. 


x,y,» 


P NU XY [clY ST Nq + Yoc + A ) 


(5.13) 


for the perturbation rotation parameters; 


1*; =fip < 

<dNU\ . , 

^dA - ) XY ( arsT ”° + Tgc + 

(5.14) 

w= np < 

<dNU\ VT , t . . 

^ dB- ) ^ Y ( arsTN » Tgc + 

(5.15) 

for the nutation series amplitudes; 




^ dAip ) ^ Y '• aTSTN ° + Tgc + 

(5.16) 


f C^) XY ( aTsT No +rcc + A) 

(5.17) 


for the nutation offsets; 

d{UTl - UTC) = ^ PN d{UTl - ufc] XY ^ aTsTNo + T ° c + ^ 5 ‘ 18 ^ 

for universal time, and 


= n PNU (H^) (pcs st No + rcc + A) (5.19) 

for polar motion. The only quantities remaining to be evaluated in Eqs. (5.12— 5.19) are the partial 
derivatives of the various rotation matrices. These are exactly the same as the partials entering 
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modeling of VLBI observations, and are set out in some detail in the publication of Sovers and 
Fanselow (1987), Sec. 2.12. For example, the first partial in Eq. (5.12) is 


4S-=fo 


Continuing the evaluation of partial derivatives with respect to geometric parameters, and using 
Q = VtPNUXY from Eq. (2.18) to simplify notation, we have for the station coordinates and their 
time rates of change, 


TSTNo = TsTN + — to) 

(5.21) 

rsTN = Q{otrsTN 0 + r GC + A) 

(5.22) 

or STN 

(5.23) 

MO 

ot STN 

(5.24) 

For the geocenter offset components, 

dr S TN q 

dr GC V 

(5.25) 

For the coordinate scale factor, 

d*STN ~ 

da = <?r5 ™° 

(5.26) 


0 0 \ 

0 -1 (5.20) 

0 0 j 


The only parameters in the tidal displacement A are the two solid-Earth-tide Love numbers h 
and /, and the tide phase angle ip. The corresponding partial derivatives of the range are somewhat 
complicated: 


dr S TN 

dh 

&tstn 

dl 

drsTN 


= QVW 


= QVW g 2$ 

\gs*. 


= QVW 


d'l’i 

dg 2 (*') 

dj> " drf>i 

dg 3 (») 

V dxjji 

where the tide phase partials of the g* s of Eqs. (2.6— 2.8) are 


dg2. 3 \t.r% 


dip 


R* 


dgi> 

dip RS 

\*p\ 


^^-hr p R.{X.y p -Y.x p ) 




(5.27) 

( 5 . 28 ) 


(5.29) 


(r p • R,)(x p X. + y p Y,) - {Y,x p - X.y p f 


(5.30) 

(5.31) 
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dg 3 . _ , , v , 

di> ~ R* ^ x,yp Y,xp ) 


J x l + vl z . , = (2r p • R, - z p Z,) 

\Ap + Vp 


(5.32) 


Finally, the partial derivative of range with respect to the post-Newtonian parameter 'yppu [see 
Eq. (3.8)] is 

— ^ ln r ‘ + r P (5.33) 

dlPPN c + r 2 — ri 2 


5.2 CLOCK PARAMETERS 

Clock parameters include the coefficients of the station and space vehicle clock error models and 
also the pseudo-range and carrier phase biases. From Eq. (5.1), the partial of pseudo-range with 
respect to the first station clock error parameter is 

_gg— = c _* — ( t3-f 3 ) + 

da>STN,R da>STNtJl dt 3 dasrN, z 

= C - p (5.34) 


where we have used 


d&STN,Z 

dt$ d 


dasTN,Z dasTN^Z 


(*3 ~t 3 ) = 1 

[^3 “ (*3 ~ ^3)] = 


and 


i> = ^~ 

* dt 3 


The range rate p is, from Eq. (5.6), given by 


(5.35) 

(5.36) 


(5.37) 


(5.38) 


The partial on the right hand side is obtained from implicit differentiation of Eq. (5.7): 
2c 2 (t 3 - t 2 ) (l - = 2 ( TSTN - Tsv) ■ (tstn - isv 


(5.39) 


dt 2 

Solving for — — gives 
dt 3 


dh __ cp — (tstn ~ rsv ) • tstn 
dt 3 cp — ( tstn — rsv) * Tsv 


where we have used Eq. (5.6). Substituting into Eq. (5.38) gives 


P = 


(fstjy — rsv) • (tstn ~ rsv ) 



tstn ~ tsv 
P 



(5.40) 


(5.41) 


The station velocity is obtained from the time derivative of Eq. (5.5), and the space vehicle velocity 
is obtained from interpolation of the PV file. 
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The partial with respect to the first space vehicle clock error parameter is 


dZ c 


— c 


a 


d<*sv,z dasv,z 


-(ti — t 2 ) ~ c 


The pseudo-range bias partial is 


dZ c 

SBstn 


= c 


( 5 . 42 ) 


( 5 . 43 ) 


The pseudo-range partials with respect to the other clock error parameters are similarly computed, 
as are the carrier phase partials; here we summarize the results. 


do.STN,R 

dbsTN,R 

dZ c 

dcsTN,R 


= C-p 

= (c — p) (t 3 ~* 3o ) 
= (c-p)(f 3 -? 3o ) 2 

dZP_ 

dBsTN 


dRP_ 

dasv,R 

d _ ZF_ 

dbsv,z 

d _ ZP_ 

dcsv,R 


= c 


= < ; (*2 — < 2 0 ) 
= c(t 2 - t 2o ) 2 


( 5 . 44 ) 


= c 


( 5 . 45 ) 


dy> c 

dasTN.v 

d y> c 
dbsTN,p 

d<p c 

dcsTN,v 


(c-p) 

(c - p)(f 3 - t 3o ) 
(c-p)(f 3 -f 3o ) 2 


chpC 

dasv, v 

dy?_ 

dbsv, v 

d<p G 

9 csv,ifi 


= c 


= C{i 2 - t 2o ) 

= c(i 2 - t 2o ) 2 


dBsv.STN 


( 5 . 46 ) 


( 5 . 47 ) 
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5.3 TROPOSPHERE PARAMETERS 


Partials of the range with respect to the dry and wet zenith delays are obtained from Eqs. (4.3) 
and (4.4): 

= [F d ry(E) + 2p Zirv F bl {E)/A]/smE 

VpZdry 

+ [2pz„ et Fb2(E)/ A + 3p 2 z dry Fb*(E)/ A 2 ]/ sin E (5.48) 

~ — - — = [F we t{E) -f 2pz dry Fb2{E)/ & + 2pz wei Fb3{E)/ A]/ sin E (5.49) 

VPZwt t 


For the zenith rate parameters: 


dp 

d PZd,v 


(t ^o)^d,w 


(5.50) 


In the analysis of data for which meteorological parameters are not available, it is convenient 
to introduce an approximation to the mapping function [Eqs. (4.3) and (4.4)] which involves a 
one-parameter estimate. This parameter p accounts for deviations from standard meteorological 
conditions. The tropospheric range is expressed as 


I 

1 

where the partial derivative is 


P = [pz dry +PSwet)/ sin ' E? + P 


dp_ 

dp 


dp (PZdry + Pz wet )uMno 

dp G(Mno i u) [l •+■ G(Mno t u)] sin E 

PZwet^^llO — Afioi/Mooi) 

G(Muo, u)G (Miox / Moot) u) [C?(Miio, u ) + (2(Mioi/A/ooi, tx)] sin E 


(5.51) 


(5.52) 
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5.4 SATELLITE PARAMETERS 


The final class of parameters in modeling GPS range measurements includes the dynamic pa- 
rameters characterizing the forces on the satellites, and the six “epoch state” position and velocity 
parameters, Tsv 0 > *sv 0 for each space vehicle. Carrier phase and pseudo-range are affected in the same 
way, through the geometric range p, by these parameters. Here we use £ to represent any parameter 
of this class. The partial derivative of geometric range with respect to £ is calculated from Eq. (5.6): 


dp dt 2 


(5.53) 


The partial on the right hand side is obtained from implicit differentiation of Eq. (5.7): 

2 c 2 ( t 3 “ h) (— - gj ) = 2(tstn ~ (~^ r 5 V r ) 

= “2(r stn ~ rsv) T (5.54) 


Solving for gives 

d£ 


dt 2 

d£ 


(rsTN - *sv) 


cp 


(*- 


* d*sv 

dtu 


*STN — Tsv 


r -f) 


Substituting this result into Eq. (5.53) gives 


dp 


— ( : t STN — r 5 y) r 


,(l- 


d*sv 

-K- 


TSTN — TSV TSV \ 
P C ) 


(5.55) 


(5.56) 


The space vehicle velocity and the variational partial 

t 2 * Detailed descriptions of the spacecraft force models and 
Mathematical Description document (Wu et al 1986). 


drsv I 

— — — are read from the PV file at time 

d£ It, 

variational partials are found in the OASIS 
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APPENDIX A 


Parameter Units and Physical Constants 


Early during the development of software for the GIPSY project, it was decided that the units 
of all parameters and partial derivatives shall be kks (kilometers, kilograms, seconds), instead of the 
more customary mks. With the exception of clock parameters and range biases, we adhere to the kks 
convention. To prevent numerical instabilities in parameter estimation arising from the disparity of 
magnitudes of the clock parameters and the observables, clock offsets and range biases are in units of 
/zsec (1CT 6 sec), rates in 10” 12 sec/sec, and rate rates in 10~ 18 sec/sec 2 . All lengths were measured 
in units of km, velocities in km/sec, and accelerations in km/sec 2 . Angular quantities, such as tide 
lag, UTPM, the perturbation tweaks, and nutation angles are in radians (or rad/sec for their rates). 
Finally, the zenith troposphere quantities are measured in km, and their rates in km/sec. These 
conventions do not necessarily apply to the formatted GPSOMC input file, where units may be used for 
some parameters that make them more easily recognizable. 

We have tried to use the constants recommended by the IAU project MERIT (Melbourne et al. f 
1983). Those that have not been defined in the text above, but which affect the results that are 
obtained using GPSOMC, are given below: 


Symbol 

Value 

Quantity 

c 

299792.458 

Velocity of light (km/sec) 

Re 

6378.140 

Equatorial radius of the Earth (km) 

OJ E 

7.2921151467 x 10“ 5 

Rotation rate of the Earth (rad/sec) 

f 

298.257 

Flattening factor of the geoid 

h 

0.609 

Vertical Love number 

l 

0.0852 

Horizontal Love number 
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APPENDIX B 


GPSOMC Parameters 


Table B.I 

Glossary of GPSOMC Parameters 


Parameter 

GPSOMC name 

Definition 

Reference 

m 

SAT PEP 0 www 

Coefficients in 

(3.20) 


SAT PRATvvwvv 

space vehicle clock 

(3.20) 


SAT PRRTvvvvvv 

model (phase) 

(3.20) 

asTN,v 

STA PEPO ss 

Coefficients in 

(3.19) 

t>STN,<p 

STA PRAT ss 

station clock 

(3.19) 

CSTN,<p 

STA PRRT ss 

model (phase) 

(3.19) 

&SV,STN 

BIAS ~vv_ss ; nn 

Carrier phase bias 

(3.18) 

<*SV, k 

SAT REPOvvvvvv 

Coefficients in 

(3.12) 

bsv, R 

SAT RRATvvvvw 

space vehicle clock 

(3.12) 

CSV,R 

SAT RRRTvvww 

model (range) 

(3.12) 

&STN, R 

STA REPO ss 

Coefficients in 

(3.4) 

t>STN,Z 

STA RRAT ss 

station clock 

(3.4) 

CSTN.Z 

STA RRRT ss 

model (range) 

(3.4) 

Bstn 

BIAS PSR ss 

Pseudo-range bias 

(3.3) 

asv 

SAT CEPOvvvvvv 

Coefficients in 

(3.22) 

bsv 

SAT CRATvvvvvv 

space vehicle 

(3.22) 

CSV 

SAT CRRTvvvvvv 

clock model 

(3.22) 

a STN 

STA CEPO ss 

Coefficients in 

(3.21) 

bsTN 

STA CRAT ss 

station clock 

(3.21) 

CSTN 

STA CRRT ss 

model 

(3.21) 


vvvvvv satellite name 
vv satellite number 

ss station number 

nn sequence number 
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Table B.I cont. 

Glossary of GPSOMC Parameters 


Parameter 

GPSOMC name 

Definition 

Reference 

r *p 

RSPINAX ss 

Cylindrical 

■RHI 

X 

LONGTUD ss 

station 


z 

POLPROJ ss 

coordinates 


'•p 

DRSP/DT 88 

Time rates of 

(2.1) 

X 

DLON/DT ss 

change of 

(2.2) ! 

i 

DPOL/DT ss 

station coordinates 

(2.3) 

h 

VLOVE ss 

Vertical Love number 

(2.5) 

i 

HLOVE ss 

Horizontal Love number 

(2.5) 

V> 

TIDPHAS ss 

Tide lag 

(2.10) 

1PPN 

GEN REL GAMMA 

Post-Newtonian gamma 

(3.8) 

XGC 

GEOCENTER X 

Coordinate frame offset 

(2.73) 

yac 

GEOCENTER Y 

Cartesian 

(2.73) 

Z GC 

GEOCENTER Z 

components 

(2.73) 

a 

COORD SCALE 

Scale factor 

(2.75) 

xsv 0 

X vvvvvv 

Space vehicle epoch 


ysv 0 

Y vvvvvv 

position Cartesian 


z SVo 

Z vvvvvv 

components 

Sec. 5.4 

isv a 

DX vvvvvv 

Space vehicle epoch 


ysvo 

DY vvvvvv 

velocity Cartesian 


ZSVo 

DZ vvvvvv 

components 



ss station number 

vvvvvv satellite name 
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Table B.I cont. 

Glossary of GPSOMC Parameters 


Parameter 

GPSOMC name 

Definition 

Reference 

©i 

X POLE MOTION 

Pole position 

(2.21) 

©2 

Y POLE MOTION 

components 

(2.20) 

UT1 - UTC 

UT1 MINUS UTC 

UTl - UTC 

(2.24) 


t ROT TWEAK OFFS 

Perturbation rotation 

(2.70) 

6®x,y,z 

t ROT TWEAK RATE 

coefficients 

(2.70) 

^oy, M y 

NUT AMPLPSI ynnn 

Nutation amplitudes 

(2.45, 2.47) 

^iy » ^3y 

NUT AMPLPSITynnn 

in longitude 

(2.45, 2.47) 

Arp 

NUT AMPLPSIA 

Longitude nutation tweak 

(2.51) 

^oy, B 23 

NUT AMPLEPS ynnn 

Nutation amplitudes 

(2.46, 2.48) 

Bij> B$j 

NUT AMPLEPS Tynnn 

in obliquity 

(2.46, 2.48) 

Ac 

NUT AMPLEPSA 

Obliquity nutation tweak 

(2.52) 

PZ&ry 

DRYZTROP as 

Dry zenith delay 

(4.1) 

PZwet 

WETZTRQP as 

Wet zenith delay 

(4.1) 

PZ&ry 

DDTRP/DT ss 

Zenith delay 

(4.2) 

pz wet 

DWTRP/DT ss 

time rates 

(4.2) 

p 

DRYMAPSG ss 

Lanyi map parameter 

(5.51) 


t X, Y, or Z 

y S or C for sine or cosine 

mm number of term in Wahr series 

sa station number 


37 









